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Abstract. Based on primitive model computer simulations with explicit microions, 
we calculate the effective interactions in a binary mixture of charged colloids with 
species A and B for different size and charge ratios. An optimal pairwise interaction 
is obtained by fitting the many-body effective forces. This interaction is close to a 
Yukawa (or Derjaguin-Landau-Verwey-Overbeek(DLVO)) pair potential but the AB 
cross-interaction is different from the geometric mean of the two direct AA and BB 
interactions. As a function of charge asymmetry, the corresponding nonadditivity 
parameter is first positive, then getting significantly negative and is getting then 
positive again. We finally show that an inclusion of nonadditivity within an optimal 
effective Yukawa model gives better predictions for the fluid pair structure than DLVO- 
theory. 
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1. Introduction 

Phase diagrams and structural correlations in binary mixtures are much richer than 
those of their one-component counterparts pj since there are additional thermodynamic 
degrees of freedom. Understanding the phase behaviour from first principles [21 |3] 
requires the knowledge of the effective interaction forces between the different species 
which is - in general - a many-body force. Even if this interaction is pairwise additive, 
the full calculation of structure and phase behaviour has only been done for selected 
cases. Among those are hard spheres [H El [6l [7], oppositely charged colloids P |9], 
two-dimensional dipolar mixtures flUl [TT] and two-dimensional Yukawa mixtures [12] • 

In this paper we consider a three-dimensional binary colloidal suspension of two 
species A and B of charged spheres ("macroions") with different charges {Zac and Zse, 
e denoting the electron charge) and diameters {a a and as) [l3l[Tll[T5]. The traditional 
Derjaguin-Landau-Verwey-Overbeek (DLVO) theory describes the interaction between 
the two species as an effective pairwise Yukawa potential [Ml [TT] 

^/ X ^ Z,e ZjC exp{-KD{r - {aj + aj)/2)) 

l + aiKD/21 + ajKn/2 er 
where [ij) = (^^4), (AB), (BB), e denotes the dielectric constant of the solvent and 
is the Debye-Hiickel screening parameter. The latter is given as 

Kn^ = ATr{J2z]pj)/ekBT (2) 
j 

where the sum runs over all microions with their charges zj and partial number densities 
Pj. The DLVO theory is a linearized theory and therefore neglects nonlinear screening 
effects [iniiTH] which give rise to effective many-body forces [20l [2TI [22l [23l [2l] . Nonlinear 
effects can at least partially be accounted for by charge renormalization which is 
conveniently calculated in a spherical Poisson-Boltzmann cell model [25]. The cell 
approach was recently generalized towards binary mixtures by Torres, Tellez and van 
Roij [26]. In the latter approach, it was shown that charge renormalization is different 
for the different species such that the ratio of effective charges is different from that of 
the bare charge. However, the cross-interaction was not addressed in this study. 

The importance of the cross-interaction between A and B relative to the direct 
part AA and BB determines the so-called non-additivity parameter A of the mixture 
which is crucial for the topology of phase diagrams. Binary hard sphere systems have 
been studied as a prototype for any non-additive mixtures [271 [2HI [29]. In general, a 
positive non-additivity is realized if the cross- interaction is more repulsive than the mean 
of the two direct interactions. For high positive non-additivity (A > 0), macrophase 
separation into an A-rich and i?-rich phase is observed, i.e. the system minimizes the 
interface where the AB cross-interactions plays a dominant role. The other case of 
negative nonadditivity (A < 0) implies a weaker cross-interaction in terms of the bare 
ones such that the system tends to mix and to exhibit micro-phase-separation [30] . 

For pairwise Yukawa interactions Vij{r) = Z*^-^ exp(— K£)r)/r {{i,j) = 
(AA), (AB), (BB)), a dimensionless nonadditivity parameter A can be quantified by 
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invoking the deviation of an ideal Berthelot mixing rule ^Tl [32] via the relation 

l + A = Z*2Vz*iZ*2 (3) 

For charged suspensions, classical DLVO theory (see eqn. (1)) implies a vanishing A 
since the effective charges are the same in all interactions. There are other realizations 
of a binary Yukawa system in dusty plasmas [331 EH [35] and metallic mixtures [36] 
or amorphous silica [37]. In fact, the binary Yukawa model has been widely used 
and employed to investigate effective interactions |38j, fluid-fluid phase separation 
[391 Uni [H], vitrification [121 SSI El HSl [16] and transport properties [17]. In most of 
studies of binary Yukawa systems, the non-additivity parameter is set to zero, except for 
Refs. [HH [32] where the effect of positive nonadditivity A on fiuid-fiuid phase separation 
is considered. In the context of dusty plasmas, there is another recent study showing 
that the non-additivity parameter A is positive in general [H]. This leads to macrophase 
separation in binary dusty plasmas, as observed in experiments [35lll8]. The physical 
origin of the interaction in dusty plasmas, however, is different from that relevant for 
charged colloidal suspensions. While for the former the ion are described by a Gurevich 
distribution, a Boltzmann distribution is appropriate for the latter. 

In this paper, we focus on the nonadditivity of the cross-interaction for charged 
colloidal suspensions. Using computer simulations with explicit microions [IHl [Ml [SH 
[52] . we calculate the effective interactions in a charged binary mixture and find that 
the sign of the non-additivity depends on the parameters, in particular on the charge 
asymmetry. The nonadditivity parameter is calculated first by using simulations of three 
pairs of macroions, namely AA, AB, and BB in a periodically repeated simulation box 
at fixed screening. We also consider larger systems with 24 macroions at different 
compositions in order to check the effect on many-body forces on A. In the latter case 
we fit the many-body forces by effective pairwise forces —dVij{r)/dr and extract A from 
the optimal fit [53]- We confirm that A is unchanged. Our main findings are i) that 
A is typically large and cannot be neglected and ii) that the sign of A depends on the 
parameters as e.g. charge asymmetry. If a binary charged colloidal mixture is described 
by an effective Yukawa model, A needs to be incorporated into the description. For 
instance when compared to DLVO theory a much better description of the fluid pair 
structure is achieved within an effective Yukawa model and non-vanishing A. 

The paper is organized as follows: In Sec. II we describe the model and the 
simulation method and apply it to the case of two macroions. Many-body simulations 
with 24 macroions are discussed in Sec. III. Then we present data for a large system with 
effective pairwise Yukawa forces in order to see the effect of non- vanishing nonadditivity 
in Sec. IV. Finally we conclude in Sec. V. 

2. Simulations with two macroions 

We model all ions as uniformly charged hard spheres such that they are interacting via 
excluded volume and Coulomb forces which are reduced by the dielectric constant e of 
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the solvent. The two species of charged colloids have a mesoscopic hardcore diameter 
aA (resp. 0"^) and a total charge Z^e (resp. Zbg) while all microions are monovalent 
with a charge e (e denoting the electron charge) and a microscopic hardcore diameter 
(Jc. For finite salt concentrations, there are both counter- and colons in the solution and 
the microscopic core of oppositely charged microions is needed to prevent the system 
from the Coulomb collapse. The averaged concentration of added salt is denoted with 
Ug. The salt is always monovalent. The system is kept at room temperature T such 
that the Bjerrum length for the microions is = e^/eksT = 7.8A with e = 80 the 
dielectric constant of water at room temperature. 

A cubic simulation box of edge length L with periodic boundary conditions is 
used containing two macroions and the following three cases are studied separatedly: i) 
two A macroions, ii) two B macroions, and iii) one A and one B macroion. The two 
macroions are placed along the room diagonal of the simulation box and possess a fixed 
central distance r. At fixed macroion positions, the microions are moved by constant 
temperature molecular dynamics and the averaged force F acting on the two macroions 
is calculated. The latter fulfills Newton's third law. For more technical details we refer 
to the Refs. |l9l [50l [SH ES, [56l [57j. Then the distance r is varied and force-distance 
curve F{r) is gained. 

Inspired by the DLVO expression (1), we anticipate that the screening length will 
not differ much in the three cases i), ii), iii) and that it will be comparable to the 
Debye-Hiickel expression. This assumption will be tested and justified later by a many 
macroion simulation reported in Section III. Therefore we adjust the box length in the 
three cases in order to reproduce the same Debye-Hiickel screening length (2). 

Results for the distance-resolved forces F{r) are shown in Figures [lH3] for the 
three cases i), ii), and iii) for extremely dilute macroion suspension with a packing 
fraction rj = 0.005. Such dilute case is chosen to diminish the boundary effects 
of the simulation box on the interaction forces. The parameters used are (Tyi=1220 
A, (Tb=680A, Za = 580, and Zb = 330. The prescribed screening is ki;)0'a= 0.8 
corresponding to box lengths of L = 6.24o"yi (for i)), L = 5.72aA (for ii)), and L = G.Octa 
(for iii)). The obtained data for the distance-resolved forces F{r) were fitted with the 
Yukawa expression Zj*-^e^/(er)(l/r + K)exp(— Kr) with (ij) = (AA), (AB), (BB). The 
screening parameter k and the three effective charge numbers Z"^^, Z^^, and Z\^ 
are used as fit parameters. We obtain kct^ = 0.81 (very close to its Debye-Hiickel 
expression) and effective charge numbers of Z\j=A7Q, Z^g=260, and Z^^=330 such 
that the nonadditivity parameter is negative: 

A = {Z\Bf /{ZUZ*bb)-1 = -^AI (4) 

In a similar manner we calculate the interaction forces and Yukawa fitting 
parameters for another packing fraction r] = 0.017 at which the images of macroions 
in neighboring cells start to affect the long-range macroion-macroion interactions. We 
obtain kdct^ = 1.15 and effective charge numbers of Z\^=505, Z|j^=265, and Z^^=342 
such that the nonadditivity parameter is A=-0.13. 
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The results for extreme dilute t] = 0.005 and dilute rj = 0.017 cases shown in 
Figures [TH3] reveal that first of all, the Yukawa expression for the forces is an excellent 
fit over the relevant distance range explored. Moreover, while the screening parameter 
is very close to its Debye-Hiickel expression, the nonadditivity of 11-13 percent is 
significant. 

Next we explore the dependence of A on the charge asymmetry a = Z^/Zb by 
changing it in the range from to 1 while keeping the size asymmetry (7^1/0"^= 1.8 
unchanged. In detail, we consider the 5-charges Zb=0, 100, 330, 580 with fixed 
Za = 580 and K£)Cr^=0.8. Using the same simulation technique and fitting procedure, we 
extract the nonadditivity parameter A according to Eq.(jl]). The results are presented in 
Figure HI The nonadditivity parameter A shows a clear non-monotonicity as a function 
of Zb starting from positive values and turning to negative ones and back to positive ones 
as Zb is increasing. For Zb = 0, the BB interaction is small, but the AB interaction has 
still a repulsive contribution from entropic contact force [51] which drives A altogether 
towards a positive value. The other cases are less intuitive and we do not have a simple 
argument for the sign of A. 

We have finally considered also the case of size-symmetric {a a = cib) but charge- 
asymmetric Za/ Zb=^-'JQ macroions at a fixed screening length of kd^a = 0.8. Here, 
the nonadditivity parameter A was found to be -0.01, much smaller than for the 
corresponding size-asymmetric case with o"a/o"b = 1.8. Hence size-asymmetry appears 
to be the more crucial input for the nonadditivity. 

3. Simulations with many macroions 

Let us now turn to a many body simulation of the primitive model in a cubic cell 
containing altogether = A^i + Nb = 24 macroions with different compositions 
X = Nb/{Na + Nb). The simulation box contains also Nc = NaZa + NbZb oppositely 
charged counterions, and Ng salt ion pairs. The other parameters are as before if not 
otherwise stated. Six different macroion packing fractions were considered: 7]= 0.017, 
0.034, 0.12, 0.16, 0.23, 0.3. For all simulations the salt concentration was kept constant 
at ?T,s=4xlO~^ mol/1. Now both the microions and the macroions are moved by constant 
temperature molecular dynamics. A simulation snapshot is shown in Figure [51 After 
equilibration, we stored 200 statistically independent macroion configurations for each 
run. For each stored configuration, we averaged the total forces Fi (i=l, A^) acting 
on the ith macroion over the microionic degrees of freedom. These forces are clearly 
many-body forces, in general. Following the idea of Ref. |53], we assign an optimal 
effective pair interaction by fitting all forces Fi in all stored configurations by the same 
pairwise Yukawa interaction. As in Section II, the four fitting parameters are kd and the 
three effective charge numbers Z\j^, ^*bbi ^^"^ ^*ab which determine the nonadditivity 
A directly. 

The results of the optimal fit for A, Z*^^, and kdCt^ are shown in Figures [6H9] 
as a function of the varied macroion volume fraction 77 for three different compositions 



Nonadditivity in binary charged colloidal suspensions 



6 



X = 1/3,1/2,2/3. 

The nonadditivity shown in Figure [6] is clearly negative and decreases with 
increasing packing fraction. This trend can be intuitively understood since asymmetries 
are amplified if one approaches smaller interparticle distances. A second important 
conclusion from Figure [6] is that the many-body simulations yield the same value for A 
as those obtained from the simulations of pairs in Sect. II. In fact, the value A = —0.13 
is reproduced at low volume fractions rj = 0.017. The effective charges and screening 
length deduced from pair macroion simulations for rj = 0.017 perfectly fit the simulation 
data for many macroions at the same packing fraction rj for the macroions. The effective 
charge numbers Z'^^^ and shown in Figures [7] and M increase slightly with volume 
fraction which is the standard trend also for one-component charged suspensions [58| l59]. 
The screening constant shown in Figure [9] increases with rj following the same trend as 
its Debye-Hiickel expression which is also indicated in Figure |9l 

4. Simulations using the optimal effective Yukawa interaction 

We finally explore the impact of a non-vanishing A on the fiuid structure of binary 
charged suspensions. In doing so we make use of the optimal effective Yukawa fit gained 
in Sect. Ill and use it as an input in classical coarse-grained binary Yukawa simulations 
(without any microions). These simulations can be done for much larger systems and we 
included 1000 A and 1000 B particles at equimolar composition. Two volume fraction 
are considered, a dilute system with 7^=0.034 and a dense system where rj=0.3. 

The effective macroion charge numbers, the screening constant and the 
nonadditivity parameter were chosen from results described in the previous section. 
In detail, for the dilute system: Z\^= 530, Z*qq= 269, fi;z)0"A=l-3, A=-0.14, while for 
the dense system, Z^^=580, Z^^=330, kz)0"a=3.4, A=-0.2. At these parameters the 
system is in the fluid phase. 

We have calculated the partial pair correlations qaa, Qbb and qab for the large 
Yukawa substitute system and the smaller system with explicit microions. The results 
are presented in Figures [10 and [Til There is good agreement for the dilute case 
demonstrating that the Yukawa flt is reproducing the pair correlations. For the dense 
system, there is again agreement except for the location of the flrst peak in 5'BB(r). This 
could have to do with flnite-size effects of the = 24 primitive model systems which 
are expected to get more prominent at high packing fractions. 

As a reference we have also performed Yukawa simulations with effective 
interactions based on two effective charges Z\j^ and Z*^^ but where A is set to zero, 
i.e. where Z\q'^ = Z^^Z^^. The differences in the pair structure are pointing to the 
importance of nonzero nonadditivity. There are even more deviation of the fluid pair 
structure when the simple DLVO expression is taken, which signiflcantly overestimates 
the structure, see the dotted lines in Figures [lO] and [Til One of the main conclusions 
therefore is that one has to be careful if A for binary charged suspensions is neglected. 
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5. Concluding remarks 

In conclusion we have determined the non-additivity parameter in binary charged 
suspensions by primitive model computer simulations with explicit microions and found 
significant deviations from zero. The sign depends on the actual parameter combination, 
in particular on the charge asymmetry. This implies that a realistic modelling of charged 
suspensions on the effective pairwise Yukawa level should incorporate a non-vanishing 
A. 

A priori an intuition for the sign of A is difficult. Intuition only works in limiting 
cases. For example, in the depletion limit of many small macroions and few big ones, 
there is attraction between the big ones which would result in a positive A. 

In the future, more detailed investigation are planned to explore the full parameter 
space better. It would be interesting to calculate the effective interaction in the solid 
phase relative to the fluid crystalline phase. 

In inhomogeneous situations like sedimentation in a gravitational fleld, binary 
suspensions have been examined by primitive model computer simulations [60] and 
a simple binary Yukawa pairwise interaction model was found to be inappropriate |61j . 
Since in our bulk simulation the pairwise Yukawa model was a good flt, we expect that 
this is due to the density gradient in the system but this will require more detailed 
studies. 

Binary suspension can be prepared in a controlled way [62l [63] and their structural 
correlation, dynamics and phase diagrams can be explored. Depending on the size and 
charge ratio different freezing diagrams (azeotropic, spindle, eutectic) |6ll EH [66], [67] 
have been obtained in the experiments. Since the nonadditivity A depends on the 
composition, this might point to the fact that details of the variation of A with 
composition and also with the phases itself (fluid or crystal) determines the shape of 
the freezing diagrams. Finally it would be interesting to study nonadditivity effects in 
charged mixtures of rods |68] or rod and sphere mixtures. It would also be interesting 
to explore the effect of multivalent counterions which could lead to mutual attraction 
[501 [69j. 
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Figure 1. (Color online) Dimensionless interaction force F\b /{kBT) for case i) 
[AA macroion pair) as a function of dimensionless separation distance t/Xb- Symbols 
denote the simulation data, the full curves are the Yukawa fit for two macroion packing 
fractions rj = 0.005 (soilid line) and rj = 0.017 (dashed line, shifted upward). The fit 
data are given in the text. 
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Figure 2. (Color online) Dimcnsionless interaction force FXb /{ksT) for case ii) {BB 
macroion pair) as a function of a dimcnsionless separation distance r/Xg. Symbols 
denote the simulation data, the full curves are the Yukawa fit for two macroion packing 
fractions rj = 0.005 (solid line) and rj = 0.017 (dashed line, shifted upward). The fit 
data are given in the text. 
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Figure 3. (Color online) Dimensionless interaction force F\b / {ksT) for case iii) 
(pair of A and B macroion) as a function of a dimensionless separation distance r/Xs- 
Symbols denote the simulation data, the full curves are the Yukawa fit for two macroion 
packing fractions -q — 0.005 (solid line) and ?/ = 0.017 (dashed line, shifted upward). 
The fit data are given in the text. 




Figure 4. (Color online) Nonadditivity parameter A for different charges Zb at a 
fixed charge Za = 580. 
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Figure 5. (Color online) Full system snapshot picture for 24 macroions with 
equimolar composition X = 1/2. The system size is 2.9 a a at a total packing fraction 
of 7^=0.3. The positively charged counter- and salt ions are shown as red dots, while 
the negatively salt ions are shown as blue dots. A vertical color gradient has been used 
for macroion positions along the 2;-axis. The parameters are: Za=580, .^^=330. 
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Figure 6. (Color online) Nonadditivity parameter A as a function of a total macroion 
packing fraction 77 for three different compositions X = 1/3, 1/2, 2/3 . 
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Figure 7. (Color online) Optimal effective AA charge number Z\j^ as a function of 
a total macroion packing fraction 77 for three different compositions X = 1/3, 1/2, 2/3. 
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Figure 8. (Color online) Optimal effective BB charge number Z'^g as a function of 
a total macroion packing fraction 77 for three different compositions X = 1/3, 1/2, 2/3. 
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Figure 9. (Color online) Optimal screening parameter kct^ as a fmiction of a total 
macroion packing fraction 77 for three different compositions X = 1/3, 1/2,2/3. The 
dashed line is the Debye-Hiickel value of screening in the simulated system according 
to Eq.(2) in the text. 
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Figure 10. (Color online) Partial macroion-macroion pair correlation functions for 
both the full primitive model (symbols) and the substitute Yukawa system (full lines) 
for a total packing fraction 77=0.034 at composition X = 1/2. The Yukawa simulations 
were carried for A=-0.14, kdo-^i^I.S, ZJj^=269. The additive Yukawa 

system results for A=0 and the same k^icta, -Z^aA' ^bb given as dashed lines. The 
DLVO predictions are included as dotted lines. 
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Figure 11. (Color online) Partial macroion-macroion pair correlation functions for 
both the full primitive model (symbols) and the substitute Yukawa system (full lines) 
for a total packing fraction 77=0.3 at composition X = 1/2. The Yukawa simulations 
were carried for A=-0.2, KDcryi=3.4, Z^^=580, Zgg~330. The additive Yukawa 
system results for A=0 and the same kd^a, Z\j^, Z^g are given as dashed lines. 
The DLVO predictions are included as dotted lines. 



